{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np  # 调用numpy\n",
    "import xarray as xr\n",
    "from eofs.standard import Eof\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def createmap():  ###############################################生成地图##########################################################\n",
    "    box = [120, 300, -20, 20]  # 经度维度\n",
    "    scale = '110m'  # 地图分辨率\n",
    "    xstep = 20  # 下面标注经纬度的步长\n",
    "    ystep = 10\n",
    "    proj = ccrs.PlateCarree(central_longitude=240)  # 确定地图投影\n",
    "    fig = plt.figure(figsize=(9, 6))  ###生成底图\n",
    "    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 确定子图，与grads的类似\n",
    "    ##海岸线\n",
    "    ax.coastlines(scale)\n",
    "    ax.set_extent(box, ccrs.PlateCarree())\n",
    "    # 标注坐标轴\n",
    "    ax.set_xticks(np.arange(box[0], box[1] + xstep, xstep), crs=ccrs.PlateCarree())\n",
    "    ax.set_yticks(np.arange(box[2], box[3] + ystep, ystep), crs=ccrs.PlateCarree())\n",
    "    # 经纬度格式，把0经度设置不加E和W\n",
    "    lon_formatter = LongitudeFormatter(zero_direction_label=False)\n",
    "    lat_formatter = LatitudeFormatter()\n",
    "    ax.xaxis.set_major_formatter(lon_formatter)\n",
    "    ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    ############################################################################################################\n",
    "    return ax, fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams['font.sans-serif'] = ['SimHei']  ###防止无法显示中文并设置黑体\n",
    "plt.rcParams['axes.unicode_minus'] = False  ###用来正常显示负号"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "iosst = xr.open_dataset('D:\\\\grads\\\\TongJi\\\\NCEP_IOSST_30y_Wt.nc')\n",
    "tpsst = xr.open_dataset('D:\\\\grads\\\\TongJi\\\\NCEP_TPSST_30y_Wt.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "tp = tpsst['st'][:, :, :]\n",
    "lat = tpsst['lat'][:]\n",
    "lon = tpsst['lon'][:]\n",
    "t = tpsst['time'][:]\n",
    "tp = tp.where(tp < 9999, np.NAN)\n",
    "avetp = tp.mean(dim=\"time\", skipna=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "jupingtp = np.zeros((len(t), len(lat), len(lon)))\n",
    "for t_i in range(len(t)):\n",
    "    jupingtp[t_i, :, :] = tp[t_i, :, :] - avetp\n",
    "a = Eof(jupingtp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "##取前几个模态\n",
    "n = 10\n",
    "eof = a.eofsAsCorrelation(neofs=n)\n",
    "pc = a.pcs(npcs=n, pcscaling=1)\n",
    "var = a.varianceFraction(neigs=n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "##\n",
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, eof[0, :, :], cmap='RdBu', transform=ccrs.PlateCarree())\n",
    "lines = ax.contour(lon, lat, eof[0, :, :], colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.2f')\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "ax.set_title('第一模态')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\eof1.1.png\")\n",
    "plt.close()\n",
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, eof[1, :, :], cmap='RdBu', transform=ccrs.PlateCarree())\n",
    "lines = ax.contour(lon, lat, eof[1, :, :], colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.2f')\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "ax.set_title('第二模态')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\eof2.1.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "##时间\n",
    "fig = plt.figure(figsize=(9, 6))\n",
    "ax = fig.subplots(1, 1)\n",
    "plt.plot(t, pc[:, 0], label='Ev=' + str(var[0] * 100)[:2] + '%')\n",
    "ax.set_yticks(np.arange(-2, 2.6, 0.2))\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "ax.set_title('PC1')\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\PC1.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(9, 6))\n",
    "ax = fig.subplots(1, 1)\n",
    "plt.plot(t, pc[:, 1], label='Ev=' + str(var[1] * 100)[:2] + '%')\n",
    "ax.set_yticks(np.arange(-2, 2.6, 0.2))\n",
    "plt.legend()\n",
    "ax.set_title('PC2')\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\PC2.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "yichangtp = tp[13, :, :] - avetp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "##合成分析\n",
    "##1983，1992 1998 为Elnino年\n",
    "##1985 1989 1999 2000 2008 为Lanina年\n",
    "tnino = np.array([4, 13, 19])\n",
    "tnina = np.array([6, 10, 20, 21, 29])\n",
    "ninos2 = np.zeros((len(lat), len(lon)))\n",
    "ninas2 = np.zeros((len(lat), len(lon)))\n",
    "avenino = np.zeros((len(lat), len(lon)))\n",
    "avenina = np.zeros((len(lat), len(lon)))\n",
    "for t_i in tnino:\n",
    "    avenino = avenino + tp[t_i, :, :]\n",
    "avenino = avenino / len(tnino)\n",
    "for t_i in tnina:\n",
    "    avenina = avenina + tp[t_i, :, :]\n",
    "avenina = avenina / len(tnina)\n",
    "for t_i in tnino:\n",
    "    ninos2 = ninos2 + (tp[t_i, :, :] - avenino) * (tp[t_i, :, :] - avenino)\n",
    "for t_i in tnina:\n",
    "    ninas2 = ninas2 + (tp[t_i, :, :] - avenina) * (tp[t_i, :, :] - avenina)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "tjianyan = (avenino - avenina) / np.sqrt(ninas2 / len(tnina) + ninos2 / len(tnino))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "###海温气候场\n",
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, avetp, cmap='bwr', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "lines = ax.contour(lon, lat, avetp, colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.1f')\n",
    "ax.set_title('海温气候场')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\qihouchang.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "###1992年El nino 海温异常场\n",
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, yichangtp, cmap='bwr', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "lines = ax.contour(lon, lat, yichangtp, colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.1f')\n",
    "ax.set_title('1992年海温异常场')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\yichang.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "### t检验\n",
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, tjianyan, cmap='bwr', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "lines = ax.contour(lon, lat, tjianyan, colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.1f')\n",
    "ax.set_title('t检验')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\tjianyan.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "##第二问1992年/1983年\n",
    "X = np.zeros((len(lat), len(lon)))\n",
    "for i in range(10):\n",
    "    X = X + eof[i, :, :] * pc[4 , i]##13-1991"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, X, cmap='bwr', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "lines = ax.contour(lon, lat, X, colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.1f')\n",
    "ax.set_title('1983年异常重构结果')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\1991yichang_chonggou1.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ax, fig = createmap()\n",
    "colorbar = ax.contourf(lon, lat, tp[4, :, :]-avetp, cmap='bwr', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(colorbar, extendrect='True', pad=0.1, fraction=0.04, shrink=1, location=\"bottom\")\n",
    "lines = ax.contour(lon, lat, tp[4, :, :]-avetp, colors='k', linestyles='-', transform=ccrs.PlateCarree())\n",
    "plt.clabel(lines, inline=True, fontsize=8, fmt='%.1f')\n",
    "ax.set_title('1983年原始场')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\1992yuanshi1.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "###主分量分析法\n",
    "iost=iosst['st'][:,:,:]\n",
    "lat1=iosst['lat'][:]\n",
    "lon1=iosst['lon'][:]\n",
    "iost=iost.where(iost<9999,np.NAN)\n",
    "aveiost=iost.mean(dim='time')\n",
    "jupingio = np.zeros((len(t), len(lat1), len(lon1)))\n",
    "for t_i in range(len(t)):\n",
    "    jupingio[t_i, :, :] = iost[t_i, :, :] - aveiost\n",
    "a=Eof(jupingio)\n",
    "iopc=a.pcs(npcs=1,pcscaling=1)\n",
    "iovar=a.varianceFraction(neigs=1)\n",
    "fig = plt.figure(figsize=(9, 6))\n",
    "ax = fig.subplots(1, 1)\n",
    "plt.plot(t, iopc[:, 0], label='Ev=' + str(iovar[0] * 100)[:2] + '%')\n",
    "ax.set_yticks(np.arange(-3.6, 1.6, 0.2))\n",
    "plt.legend()\n",
    "ax.set_title('iopc1')\n",
    "plt.savefig(\"D:\\\\grads\\\\TongJi\\\\ex6\\\\iopc1.png\")\n",
    "plt.grid()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
